#include "lomb/period.h"

#include <stdio.h>
#include <stdlib.h>

/******************************
Usage: lomb MIN-FREQ MAX-FREQ OVERSAMPLING-RATIO INPUT-FILE > OUT-FILE

OVERSAMPLING-RATIO determines the number of output freqs, as
   #of frequencies=0.5*ofac*(hifac - lowfac)*N
   where hifac = f_hi/f_nyquist and N = # of input points

INPUTFILE should be in two column form: time   amplitude

OUTFILE is in the form

frequency_index  freq   amplitude   z50  z95  z99

where z95 is the z(or y)-height of the 95% confidence level etc
- to plot the significance levels on top of the plot, use
replot "outputfile" using 2:4
replot "outputfile" using 2:5
replot "outputfile" using 2:6
in gnuplot

 ******************************/



double * time = 0;   // time series array
double * amp = 0;    // array of time series amplitudes
double * px = 0;     // frequency array
double * py = 0;     // power spectrum

// variable used as to dynamically create memory to store the data (initial
// buffer size - must be a power of 2). This will increase as necessary by
// repeated doubling, depending on the length of input time series.
int nmax = 8192 * 4L;

float ave, var;   // some variables used in normalizing the spectrum ?!

char * prog_name(); // returns the program name
void help();        // prints help text
void error();       // error text
void print_title();
int read_data(FILE* file);


int main(int argc, char ** argv)
{
    /* file pointer for input data */
    FILE* file = NULL;

    int nlength;
    int N;
    int nout;
    int jmax;
    float prob;
    int maxout;
    int norm_factor;
    double f_hi;
    double f_low;
    double firsttime;
    double M, z50, z95, z99; //variables for calculating significance levels
    int nstart;

    double oversample_factor = 4;

    for (int i = 1; i < argc; ++i)
    {
        if (*argv[i] == '-')
        {
            switch (*(argv[i]+1))
            {
                case 'h':  /* help and exit */
                    help();
                    exit(0);
                    break;
                case 'x':
                    break;
                default:
                    fprintf(stderr, "%s: unrecognised option %s ignored\n",
                            "lomb", argv[i]);
                    break;
            }
        }
        else if (i == argc-4)   /* MIN FREQUENCY */
        {
            f_low = atof(argv[i]);
            fprintf(stderr,"\n\nMin freq: %f\n", f_low);
        }
        else if (i == argc-3)   /*MAXIMUM FREQUENCY*/
        {
            f_hi = atof(argv[i]);
            printf("Max freq: %f\n", f_hi);
        }
        else if (i == argc-2)  /*OVERSAMPLING FACTOR ofac*/
        {
            oversample_factor = atof(argv[i]);
            printf("Oversampling factor: %f\n", oversample_factor);
        }
        else if (i == argc-1)  /* input file name */
        {
            if ( (file = fopen(argv[i],"rt")) == NULL) {
                fprintf("%s: can't open file '%s'\n", "lomb", argv[i]);
                exit(2);
            }
        }
    }

    if (!file)
    {
        help();
        exit(1);
    }

    // === reads N data points into time[N], amp[N] arrays
    N = read_data(file);

    // === Move time origin to zero : time[0] = 0
    firsttime = time[0];
    fprintf(stderr,"Time origin: %.12g\n",firsttime);
    for (int i = 0; i < N; ++i)
    {
        *(time + i) = *(time + i) - firsttime;
    }
    fprintf(stderr,"Time Range: %.12g\n",time[N-1]);

    // calculate periodogram, returns M = effm
    M = period(time, amp, N, oversample_factor, f_low, f_hi, px, py, 64 * nmax,
            &nout, &nstart,& jmax, &prob);

    // print output
    avevar(amp, N, &ave, &var);
    maxout = nout;
    fprintf(stderr,"\n\nnout: %d\tvar: %f nstart %d\n",nout,var,nstart);
    nlength = nout - nstart;
    norm_factor=(nlength/var);

    // Calculate z values of 50%,95% and 99% certainty levels.
    z50 = log(M/0.5);
    z95 = log(M/0.05);
    z99 = log(M/0.01);

    printf("! freq amp 50 %% 95 %% 99 %% period\n");
    for (int freq = nstart; freq < maxout; ++freq)
    {
        printf("%10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n", px[freq],
                py[freq], z50, z95, z99, 1.0/px[freq]);
    }

    fprintf(stderr,"Significance probability of highest peak: Prob: %f\n",(1-prob));

    return 0;
}









//=============================================================================
//=============================================================================
//=============================================================================


char * prog_name(char *s)
{
    char *p = s + strlen(s);
    while (p >= s && *p != '/')
        p--;
    return (p+1);
}



void error(char * s)
{
    fprintf(stderr,"%s: %s\n", "lomb", s);
    exit(1);
}

static char * help_strings[] = {
        "\n usage: %s [options ...] MIN-FREQ MAX-FREQ OVERSAMPLING-RATIO INPUT-FILE\n\n",
        "  INPUT_FILE: the name of a text file ",
        "    containing a sampled time series, presented as two",
        "    columns of numbers - sample times, sample values and error\n",
        "  OPTIONS: may be any of:",
        "    -h   Print on-line help",
        "    -x   left for as template option\n",
        0
};

void help()
{
    int i;

    fprintf(stderr, help_strings[0], "lomb");

    for (i = 1; help_strings[i] != 0; i++)
    {
        fprintf(stderr, "%s\n",help_strings[i]);
    }
}


int read_data(FILE* file)
{
    unsigned long npts = 0L;
    double err;

    if ((time = (double *)malloc(nmax * sizeof(double))) == 0 ||
            (amp = (double *)malloc(nmax * sizeof(double))) == 0 ||
            (px = (double *)malloc(nmax * sizeof(double))) == 0 ||
            (py = (double *)malloc(nmax * sizeof(double))) == 0)
    {
        if (time) (void)free(time);
        fclose(file);
        error("insufficient memory");
    }

    while (fscanf(file, "%lf%lf", &time[npts], &amp[npts]) == 2)
    {
        if (++npts >= nmax)
        {
            double *xt, *yt, *w1t, *w2t;
            unsigned int nmaxt = nmax << 1;

            if (64 * nmaxt * sizeof(double) < nmax) {
                fprintf(stderr,"%s insufficient memory, tructating at row %ld\n",
                        "lomb", npts);
                break;
            }
            if ((xt = (double *)realloc(time,nmaxt*sizeof(double))) ==  NULL) {
                fprintf(stderr,"%s insufficient memory, tructating at row %ld\n",
                        "lomb", npts);
                break;
            }
            time = xt;
            if ((yt = (double *)realloc(amp,nmaxt*sizeof(double))) ==  NULL) {
                fprintf(stderr,"%s insufficient memory, tructating at row %ld\n",
                        "lomb", npts);
                break;
            }
            amp = yt;
            if ((w1t = (double *)realloc(px,64*nmaxt*sizeof(double))) == NULL){
                fprintf(stderr,
                        "%s: insufficient memory, truncating input at row %ld\n",
                        "lomb", npts);
                break;
            }
            px = w1t;
            if ((w2t = (double *)realloc(py,64*nmaxt*sizeof(double))) == NULL){
                fprintf(stderr,
                        "%s: insufficient memory, truncating input at row %ld\n",
                        "lomb", npts);
                break;
            }
            py = w2t;
            nmax = nmaxt;
        }
    }

    fclose(file);
    if (npts < 1) error("no data read");
    return (npts);
}




